logo_datactivist           logo_pleiade


Objectif : Exploration des données et modélisations

Date début de l’analyse : 17 mai 2022

Date de la dernière modification : 23 mai 2022

# Packages
packages = c("tidyverse", "glue", "kableExtra", "knitr", "plotly", "hrbrthemes", "gridExtra", "stats", "arrow", "lme4", "tidymodels", "broom.mixed", "multilevelmod", "dotwhisker", "rAmCharts")

## Installation des packages si besoin et chargement des librairies
package.check <- lapply(packages,
  FUN = function(x) {
    if (!require(x, character.only = TRUE)) {
      install.packages(x, dependencies = TRUE)
      library(x, character.only = TRUE)}})

# Import BSO complet (2013-2020)
data <- read_csv("../data/4.process/managing-variables_2nd-part/whole_bso.csv", col_types = cols(author_position = "c"), na = c("", "NA"))

# Création de nouvelles variables : couleur de journal finale et APC fiables
data <- data %>% mutate(oa_color_finale = case_when(oa_color_article_BSO == "diamond" ~ "diamond",
                                                    oa_color_openalex == "gold" ~ "gold",
                                                    oa_color_openalex == "hybrid" ~ "hybridOA",
                                                    TRUE ~ NA_character_), #qs Maurits : catégorie "other" ou tout en NA ? Sur quelle variable se baser alors ?
                        amount_apc_EUR_trusted = case_when(apc_source == "openAPC_estimation_publisher" | apc_source == "openAPC_estimation_publisher_year" ~ NA_real_,
                                                           TRUE ~ amount_apc_EUR))

I - Analyse exploratoire


Y moyen par année et par variable catégorielle


Y1 = amount_apc_EUR

Discipline

table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR)) %>%  group_by(year, bso_classification) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR, na.rm = T)) %>% ungroup()

  # Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = bso_classification, colour = bso_classification)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  #scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par discipline et par année") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph)


Tier

table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR)) %>% group_by(year, tier) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR, na.rm = T)) %>% ungroup() %>% mutate(tier = as.character(tier))

  # Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = tier, colour = tier)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  #scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par tier et par année") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph)


French CA

table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR)) %>% 
    group_by(year, is_french_CA_bso_wos_openalex_single_lang) %>% 
    summarise(montant_moyen = mean(amount_apc_EUR, na.rm = T)) %>% ungroup() %>% 
    mutate(is_french_CA = as.character(is_french_CA_bso_wos_openalex_single_lang))

  # Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = is_french_CA, colour = is_french_CA)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  #scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si le CA est français") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph)


Couleur OA

table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR)) %>% group_by(year, oa_color_finale) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR, na.rm = T)) %>% ungroup() %>% mutate(oa_color_finale = as.character(oa_color_finale))

  # Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = oa_color_finale, colour = oa_color_finale)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  #scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par couleur OA et par année") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph)


Covered by Couperin

table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR)) %>% group_by(year, is_covered_by_couperin) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR, na.rm = T)) %>% ungroup() %>% mutate(is_covered_by_couperin = as.character(is_covered_by_couperin))

  # Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = is_covered_by_couperin, colour = is_covered_by_couperin)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  #scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si l'article est couvert par les accords Couperin") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph)


Y2 = amount_apc_EUR_trusted

Discipline

table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>%  group_by(year, bso_classification) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup()

  # Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = bso_classification, colour = bso_classification)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  #scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par discipline et par année") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph)


Au vu des graphiques présentés ci-dessus, nous décidons de ne modéliser que les articles pour lesquels le montant des APC est fiable, c’est-à-dire provenant de l’une des sources suivantes :

  • doaj ;
  • OpenAPC ;
  • openAPC_estimation_issn_year ;
  • openAPC_estimation_issn.

Parmi les 6 sources d’APC recensées par le BSO, nous écartons ainsi “openAPC_estimation_publisher” et “openAPC_estimation_publisher_year”.

# Sélection de variables pour les modèles
bso_mvp <- data %>% 
    select(year, bso_classification, apc_has_been_paid, tier, oa_color_finale, is_french_CA_bso_wos_openalex_single_lang, amount_apc_EUR_trusted, is_covered_by_couperin) %>% 
    mutate(bso_classification = as.factor(bso_classification), apc_has_been_paid = as.logical(apc_has_been_paid), tier = as.factor(tier), oa_color_finale = as.factor(oa_color_finale), is_french_CA = as.logical(is_french_CA_bso_wos_openalex_single_lang), is_covered_by_couperin = as.logical(is_covered_by_couperin))

# Préparation du jeu de données prêt à modéliser
bso_pop <- bso_mvp %>% 
    filter(!is.na(bso_classification),
           !is.na(oa_color_finale),
           !is.na(is_french_CA_bso_wos_openalex_single_lang),
           !is.na(amount_apc_EUR_trusted),
           !is.na(is_covered_by_couperin)) %>% 
    filter(apc_has_been_paid == TRUE, amount_apc_EUR_trusted > 0)

À partir de cela, plusieurs manipulations sont réalisées pour préparer les données aux modélisations :

  • sélection : sélection des variables qui entreront dans les modélisations ;
  • filtres : premier filtre pour ne garder que les articles pour lesquels toutes les variables à intégrer aux modèles ont une valeur connue (not NA), puis second filtre sur les articles pour lesquels des APC ont été payés (apc_has_been_paid == TRUE, amount_apc_EUR_trusted > 0).

Le nombre d’articles dans notre base finale pour les modélisations est ainsi de 119616. En ne gardant que les sources d’APC fiables, 28554 articles ont été écartés, soit 19.27% des articles avec des APC payés et valeurs connues pour les autres variables.


Tier

table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>% group_by(year, tier) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>% mutate(tier = as.character(tier))

  # Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = tier, colour = tier)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  #scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par tier et par année") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph)


Au vu des graphiques présentés ci-dessus, nous décidons de ne modéliser que les articles pour lesquels le montant des APC est fiable, c’est-à-dire provenant de l’une des sources suivantes :

  • doaj ;
  • OpenAPC ;
  • openAPC_estimation_issn_year ;
  • openAPC_estimation_issn.

Parmi les 6 sources d’APC recensées par le BSO, nous écartons ainsi “openAPC_estimation_publisher” et “openAPC_estimation_publisher_year”.

# Sélection de variables pour les modèles
bso_mvp <- data %>% 
    select(year, bso_classification, apc_has_been_paid, tier, oa_color_finale, is_french_CA_bso_wos_openalex_single_lang, amount_apc_EUR_trusted, is_covered_by_couperin) %>% 
    mutate(bso_classification = as.factor(bso_classification), apc_has_been_paid = as.logical(apc_has_been_paid), tier = as.factor(tier), oa_color_finale = as.factor(oa_color_finale), is_french_CA = as.logical(is_french_CA_bso_wos_openalex_single_lang), is_covered_by_couperin = as.logical(is_covered_by_couperin))

# Préparation du jeu de données prêt à modéliser
bso_pop <- bso_mvp %>% 
    filter(!is.na(bso_classification),
           !is.na(oa_color_finale),
           !is.na(is_french_CA_bso_wos_openalex_single_lang),
           !is.na(amount_apc_EUR_trusted),
           !is.na(is_covered_by_couperin)) %>% 
    filter(apc_has_been_paid == TRUE, amount_apc_EUR_trusted > 0)

À partir de cela, plusieurs manipulations sont réalisées pour préparer les données aux modélisations :

  • sélection : sélection des variables qui entreront dans les modélisations ;
  • filtres : premier filtre pour ne garder que les articles pour lesquels toutes les variables à intégrer aux modèles ont une valeur connue (not NA), puis second filtre sur les articles pour lesquels des APC ont été payés (apc_has_been_paid == TRUE, amount_apc_EUR_trusted > 0).

Le nombre d’articles dans notre base finale pour les modélisations est ainsi de 119616. En ne gardant que les sources d’APC fiables, 28554 articles ont été écartés, soit 19.27% des articles avec des APC payés et valeurs connues pour les autres variables.


French CA

table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>% 
    group_by(year, is_french_CA_bso_wos_openalex_single_lang) %>% 
    summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>% 
    mutate(is_french_CA = as.character(is_french_CA_bso_wos_openalex_single_lang))

  # Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = is_french_CA, colour = is_french_CA)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  #scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si le CA est français") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph)


Au vu des graphiques présentés ci-dessus, nous décidons de ne modéliser que les articles pour lesquels le montant des APC est fiable, c’est-à-dire provenant de l’une des sources suivantes :

  • doaj ;
  • OpenAPC ;
  • openAPC_estimation_issn_year ;
  • openAPC_estimation_issn.

Parmi les 6 sources d’APC recensées par le BSO, nous écartons ainsi “openAPC_estimation_publisher” et “openAPC_estimation_publisher_year”.

# Sélection de variables pour les modèles
bso_mvp <- data %>% 
    select(year, bso_classification, apc_has_been_paid, tier, oa_color_finale, is_french_CA_bso_wos_openalex_single_lang, amount_apc_EUR_trusted, is_covered_by_couperin) %>% 
    mutate(bso_classification = as.factor(bso_classification), apc_has_been_paid = as.logical(apc_has_been_paid), tier = as.factor(tier), oa_color_finale = as.factor(oa_color_finale), is_french_CA = as.logical(is_french_CA_bso_wos_openalex_single_lang), is_covered_by_couperin = as.logical(is_covered_by_couperin))

# Préparation du jeu de données prêt à modéliser
bso_pop <- bso_mvp %>% 
    filter(!is.na(bso_classification),
           !is.na(oa_color_finale),
           !is.na(is_french_CA_bso_wos_openalex_single_lang),
           !is.na(amount_apc_EUR_trusted),
           !is.na(is_covered_by_couperin)) %>% 
    filter(apc_has_been_paid == TRUE, amount_apc_EUR_trusted > 0)

À partir de cela, plusieurs manipulations sont réalisées pour préparer les données aux modélisations :

  • sélection : sélection des variables qui entreront dans les modélisations ;
  • filtres : premier filtre pour ne garder que les articles pour lesquels toutes les variables à intégrer aux modèles ont une valeur connue (not NA), puis second filtre sur les articles pour lesquels des APC ont été payés (apc_has_been_paid == TRUE, amount_apc_EUR_trusted > 0).

Le nombre d’articles dans notre base finale pour les modélisations est ainsi de 119616. En ne gardant que les sources d’APC fiables, 28554 articles ont été écartés, soit 19.27% des articles avec des APC payés et valeurs connues pour les autres variables.


Couleur OA

table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>% group_by(year, oa_color_finale) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>% mutate(oa_color_finale = as.character(oa_color_finale))

  # Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = oa_color_finale, colour = oa_color_finale)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  #scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par couleur OA et par année") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph)


Au vu des graphiques présentés ci-dessus, nous décidons de ne modéliser que les articles pour lesquels le montant des APC est fiable, c’est-à-dire provenant de l’une des sources suivantes :

  • doaj ;
  • OpenAPC ;
  • openAPC_estimation_issn_year ;
  • openAPC_estimation_issn.

Parmi les 6 sources d’APC recensées par le BSO, nous écartons ainsi “openAPC_estimation_publisher” et “openAPC_estimation_publisher_year”.

# Sélection de variables pour les modèles
bso_mvp <- data %>% 
    select(year, bso_classification, apc_has_been_paid, tier, oa_color_finale, is_french_CA_bso_wos_openalex_single_lang, amount_apc_EUR_trusted, is_covered_by_couperin) %>% 
    mutate(bso_classification = as.factor(bso_classification), apc_has_been_paid = as.logical(apc_has_been_paid), tier = as.factor(tier), oa_color_finale = as.factor(oa_color_finale), is_french_CA_bso_wos_openalex_single_lang = as.logical(is_french_CA_bso_wos_openalex_single_lang), is_covered_by_couperin = as.logical(is_covered_by_couperin))

# Préparation du jeu de données prêt à modéliser
bso_pop <- bso_mvp %>% 
    filter(!is.na(bso_classification),
           !is.na(oa_color_finale),
           !is.na(is_french_CA_bso_wos_openalex_single_lang),
           !is.na(amount_apc_EUR_trusted),
           !is.na(is_covered_by_couperin)) %>% 
    filter(apc_has_been_paid == TRUE, amount_apc_EUR_trusted > 0)

À partir de cela, plusieurs manipulations sont réalisées pour préparer les données aux modélisations :

  • sélection : sélection des variables qui entreront dans les modélisations ;
  • filtres : premier filtre pour ne garder que les articles pour lesquels toutes les variables à intégrer aux modèles ont une valeur connue (not NA), puis second filtre sur les articles pour lesquels des APC ont été payés (apc_has_been_paid == TRUE, amount_apc_EUR_trusted > 0).

Le nombre d’articles dans notre base finale pour les modélisations est ainsi de 119616. En ne gardant que les sources d’APC fiables, 28554 articles ont été écartés, soit 19.27% des articles avec des APC payés et valeurs connues pour les autres variables.


Covered by Couperin

table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>% group_by(year, is_covered_by_couperin) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>% mutate(is_covered_by_couperin = as.character(is_covered_by_couperin))

  # Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = is_covered_by_couperin, colour = is_covered_by_couperin)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  #scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si l'article est couvert par les accords Couperin") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph)


Au vu des graphiques présentés ci-dessus, nous décidons de ne modéliser que les articles pour lesquels le montant des APC est fiable, c’est-à-dire provenant de l’une des sources suivantes :

  • doaj ;
  • OpenAPC ;
  • openAPC_estimation_issn_year ;
  • openAPC_estimation_issn.

Parmi les 6 sources d’APC recensées par le BSO, nous écartons ainsi “openAPC_estimation_publisher” et “openAPC_estimation_publisher_year”.

# Sélection de variables pour les modèles
bso_mvp <- data %>% 
    select(year, bso_classification, apc_has_been_paid, tier, oa_color_finale, is_french_CA_bso_wos_openalex_single_lang, amount_apc_EUR_trusted, is_covered_by_couperin) %>% 
    mutate(bso_classification = as.factor(bso_classification), apc_has_been_paid = as.logical(apc_has_been_paid), tier = as.factor(tier), oa_color_finale = as.factor(oa_color_finale), is_french_CA_bso_wos_openalex_single_lang = as.logical(is_french_CA_bso_wos_openalex_single_lang), is_covered_by_couperin = as.logical(is_covered_by_couperin))

# Préparation du jeu de données prêt à modéliser
bso_pop <- bso_mvp %>% 
    filter(!is.na(bso_classification),
           !is.na(oa_color_finale),
           !is.na(is_french_CA_bso_wos_openalex_single_lang),
           !is.na(amount_apc_EUR_trusted),
           !is.na(is_covered_by_couperin)) %>% 
    filter(apc_has_been_paid == TRUE, amount_apc_EUR_trusted > 0)

À partir de cela, plusieurs manipulations sont réalisées pour préparer les données aux modélisations :

  • sélection : sélection des variables qui entreront dans les modélisations ;
  • filtres : premier filtre pour ne garder que les articles pour lesquels toutes les variables à intégrer aux modèles ont une valeur connue (not NA), puis second filtre sur les articles pour lesquels des APC ont été payés (apc_has_been_paid == TRUE, amount_apc_EUR_trusted > 0).

Le nombre d’articles dans notre base finale pour les modélisations est ainsi de 119616. En ne gardant que les sources d’APC fiables, 28554 articles ont été écartés, soit 19.27% des articles avec des APC payés et valeurs connues pour les autres variables.


Valeurs aberrantes

graph <- bso_pop %>% ggplot(aes(x = amount_apc_EUR_trusted)) +
      geom_boxplot(fill = "#112446") +
      coord_flip() +
      labs(x = "Montant des APC fiables", y = " ", title = "Boxplot de Y") +
      theme_ipsum() +
      theme(axis.title.y = element_text(size=12), axis.title.x = element_text(size=12))
#boxplot <- ggplotly(graph)
graph1 <- bso_pop %>% ggplot(aes(x = amount_apc_EUR_trusted)) +
      geom_histogram(bins = 30L, fill = "#112446") +
      labs(x = "Montant des APC fiables", y = "Nombre d'articles", title = "Histogramme de Y") +
      theme_ipsum() +
      theme(axis.title.y = element_text(size=12), axis.title.x = element_text(size=12)) +
      stat_function(fun = dnorm, args = list(mean = mean(bso_pop$amount_apc_EUR_trusted), sd = sd(bso_pop$amount_apc_EUR_trusted)), col = "white")
#histo <- ggplotly(graph1)
grid.arrange(graph, graph1, ncol = 2, nrow = 1)

L’étendue (écart entre les valeurs minimales et maximales) du montant des APC est importante, avec des valeurs allant de 2.27 à 9848.4 euros. Nous constatons 2 valeurs extrêmes sur le boxplot, où les montants sont respectivement de 9848.4 et 9028.43 euros, nous filtrons donc la base sur les montants inférieurs à 7000 euros (le 3e APC le plus élevé est de 6848.67 euros) et regardons à nouveau les graphiques d’évolution des APC par année, selon les différentes variables qualitatives.


Graphiques d’évolution

Discipline

# Suppression des valeurs aberrantes
bso_pop <- bso_pop %>% filter(amount_apc_EUR_trusted < 7000) %>% mutate(is_french_CA_bso_wos_openalex_single_lang = as.logical(is_french_CA_bso_wos_openalex_single_lang))

# Graphique / discipline
    # table
table <- bso_pop %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>%  group_by(year, bso_classification) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup()
    # plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = bso_classification, colour = bso_classification)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  #scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par discipline et par année") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph)


Tier

# Graphique / tier
    # table
table <- bso_pop %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>%  group_by(year, tier) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup()
    # plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = tier, colour = tier)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  #scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par tier et par année") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph)


French CA

# Graphique / french CA
    # table
table <- bso_pop %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>%  group_by(year, is_french_CA_bso_wos_openalex_single_lang) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>% 
  rename(is_french_CA = is_french_CA_bso_wos_openalex_single_lang)
    # plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = is_french_CA, colour = is_french_CA)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  #scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si le CA est français") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph)


Couleur OA

# Graphique / oa_color_finale
    # table
table <- bso_pop %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>%  group_by(year, oa_color_finale) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup()
    # plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = oa_color_finale, colour = oa_color_finale)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  #scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par couleur OA et par année") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph)


Covered by Couperin

table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>% group_by(year, is_covered_by_couperin) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>% mutate(is_covered_by_couperin = as.character(is_covered_by_couperin))

  # Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = is_covered_by_couperin, colour = is_covered_by_couperin)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  #scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si l'article est couvert par les accords Couperin") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph)


Normalité de distribution

# Normalité de Y
mean <- mean(bso_pop$amount_apc_EUR_trusted, na.rm = TRUE)   
sd <- sd(bso_pop$amount_apc_EUR_trusted, na.rm = TRUE)  
ggplot(bso_pop, aes(bso_pop$amount_apc_EUR_trusted)) +
      geom_histogram(aes(y=..density..),  color="black", fill = "steelblue", alpha = 0.2) +
    geom_density( color="red", size = 1, adjust = 5) +
    stat_function(fun = dnorm, colour = "Black", size = 1, args = list(mean = mean, sd = sd))  +
    xlim(c(min(bso_pop$amount_apc_EUR_trusted)-1, max(bso_pop$amount_apc_EUR_trusted)+1))+ ggtitle("Distribution du montant des APC fiables") +
  xlab("Montant des APC fiables") + ylab("Densité") +
     theme_classic()

# Test statistique de normalité
library(DescTools)
JarqueBeraTest(ts(bso_pop$amount_apc_EUR_trusted), robust = FALSE, method = "chisq")  # Y non normal (on rejette H0)
Jarque Bera Test

data: ts(bso_pop$amount_apc_EUR_trusted) X-squared = 47076, df = 2, p-value < 2.2e-16

Le montant des APC fiables n’est pas normalement distribué.


Corrélations entre les variables


II - Modélisations Y : 2013-2020


Péparation des données aux modélisations

Pour les modélisations, nous cherchons donc à expliquer et prédire le montant des APC fiables grâce aux variables suivantes :

  • variables explicatives (indépendantes) : oa_color_finale, tier, bso_classification, is_french_CA_bso_wos_openalex_single_lang ;
  • sources de variations (sous-populations) : year, is_covered_by_couperin.
# Train test split
train_size <- floor(0.8 * nrow(bso_pop)) # 80% of the sample size

## set the seed to make your partition reproducible
set.seed(123)
train_ind <- sample(seq_len(nrow(bso_pop)), size = train_size)

train <- bso_pop[train_ind, ]
test <- bso_pop[-train_ind, ]

Ensuite, nous séparons notre base de données en 2 échantillons : un échantillon d’apprentissage (train) et un échantillon de test (test). Le premier contient 80% des données, soit 95691 articles et le second contient les 20% derniers, soit 23923 articles. Tous les modèles seront construits sur l’échantillon d’apprentissage, puis seront validés / testés sur l’échantillon test ; s’agissant de données inconnues (les modèles n’auront pas été entraînés dessus), les appliquer à cet échantillon permettra de tester leur robustesse et leur capacité prédictive.


Régression linéaire : effets fixes


  • Modèle avec seulement l’année comme variable explicative
# Reg linéaire
lm1 = lm(amount_apc_EUR_trusted ~ year, data = train)
summary <- tidy(lm1)
kable(summary, format = "html") %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), 
                full_width = F, 
                html_font = "sans-serif")
term estimate std.error statistic p.value
(Intercept) -77428.73049 2396.139753 -32.31395 0
year 39.27705 1.187785 33.06747 0

Chaque année, le montant des APC (fiables) augmente de 39 euros en moyenne.


  • Modèle avec toutes les variables explicatives
# Reg linéaire
lm2 = lm(amount_apc_EUR_trusted ~ year + tier + oa_color_finale + is_covered_by_couperin + is_french_CA_bso_wos_openalex_single_lang, data = train)
summary <- tidy(lm2)
kable(summary, format = "html") %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), 
                full_width = F, 
                html_font = "sans-serif")
term estimate std.error statistic p.value
(Intercept) -60191.26359 2256.790367 -26.67118 0
year 30.85978 1.118668 27.58619 0
tier2 -311.08861 9.693970 -32.09094 0
tier3 -440.85370 9.574222 -46.04590 0
tier4 -446.89060 10.263922 -43.53995 0
oa_color_finalehybridOA 801.69500 8.181967 97.98316 0
is_covered_by_couperinTRUE -96.96469 9.262517 -10.46850 0
is_french_CA_bso_wos_openalex_single_langTRUE -67.89926 4.874324 -13.92998 0


Modèles mutli-niveaux : effets fixes + aléatoires


Modèles de Joël

  • Premier modèle : random intercepts par discipline

La discipline est ici la source de variations du modèle, chaque discipline dispose de son ordonnée à l’origine.

# Modèle linéaire en utilisant le package lme4
model1_spec <- linear_reg() %>% 
    set_engine("lmer") 

# une ordonnée à l'origine qui peut varier d'une discipline à l'autre
model1_fit <- model1_spec %>% 
    fit(amount_apc_EUR_trusted ~ (1 | bso_classification), data = train) 

# vue générale du modèle
summary <- tidy(model1_fit)
kable(summary, format = "html") %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), 
                full_width = F, 
                html_font = "sans-serif")
effect group term estimate std.error statistic
fixed NA (Intercept) 1495.2871 105.4812 14.17586
ran_pars bso_classification sd__(Intercept) 332.7945 NA NA
ran_pars Residual sd__Observation 810.2233 NA NA 

# les variations des APC par discipline
tidy(model1_fit, effects = "ran_vals") %>%  # on calcule les valeurs des effets aléatoires
    mutate(estimate = 1777 + estimate,  # pour rajouter l'ordonnée à l'origine pour l'ensemble de la population
           ymin = estimate - 2 * `std.error`, # pour calculer des intervalles à 95%
           ymax = estimate + 2 * `std.error`) %>% 
    arrange(estimate) %>% 
    mutate(level = as_factor(level)) %>% 
    ggplot(aes(x = level, y = estimate)) +
    geom_pointrange(aes(ymin = ymin, ymax = ymax)) +
    geom_hline(yintercept = 1777, linetype = 2) +
    coord_flip()


  • Deuxième modèle : random intercepts par discipline et trend dépend de l’année

Par rapport au premier modèle, on ajoute l’année comme variable explicative du modèle - qui impacte donc l’angle de la pente de régression (mais celle-ci est la même par cluster, c’est-à-dire par discipline).

## Tentons un 2e modèle, dans lequel on suppose que le montant des APC augmente au cours du temps 

# on commence par transformer la variable année
train <- train %>% 
    mutate(year = year - 2013)

model2_fit <- model1_spec %>% 
    fit(amount_apc_EUR_trusted ~ year + (1 | bso_classification), data = train) # l'effet de l'année sur l'APC est supposé le même quel que soit la discipline

# vue générale du modèle
summary <- tidy(model2_fit) # l'APC moyen augmente de 34,6€ par an
kable(summary, format = "html") %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), 
                full_width = F, 
                html_font = "sans-serif")
effect group term estimate std.error statistic
fixed NA (Intercept) 1311.20846 103.378135 12.68361
fixed NA year 41.67462 1.174069 35.49587
ran_pars bso_classification sd__(Intercept) 325.72589 NA NA
ran_pars Residual sd__Observation 804.94652 NA NA 

# les variations des APC par discipline - ici les valeurs sont calculées pour 2013 (year = 0)
tidy(model2_fit, effects = "ran_vals") %>%  # on calcule les valeurs des effets aléatoires
    mutate(estimate = 1565 + estimate,  # pour rajouter l'ordonnée à l'origine pour l'ensemble de la population
           ymin = estimate - 2 * `std.error`, # pour calculer des intervalles à 95%
           ymax = estimate + 2 * `std.error`) %>% 
    arrange(estimate) %>% 
    mutate(level = as_factor(level)) %>% 
    ggplot(aes(x = level, y = estimate)) +
    geom_pointrange(aes(ymin = ymin, ymax = ymax)) +
    geom_hline(yintercept = 1565, linetype = 2) +
    coord_flip()


  • Troisième modèle : random intercepts par discipline et trend dépend de l’année par discipline

Par rapport au deuxième modèle, on autorise la pente liée à l’année, à varier selon la discipline (more random effects).


Modèles pour l’analyse

  • Premier modèle : random intercepts par année

L’année est ici la source de variations du modèle, chaque année dispose de son ordonnée à l’origine.

# Modèle linéaire en utilisant le package lme4
model1_spec <- linear_reg() %>% 
    set_engine("lmer") 

# une ordonnée à l'origine qui peut varier d'une année à l'autre
model1_fit <- model1_spec %>% 
    fit(amount_apc_EUR_trusted ~ 1 + (1 | year), data = train) 

# vue générale du modèle
summary <- tidy(model1_fit)
kable(summary, format = "html") %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), 
                full_width = F, 
                html_font = "sans-serif")
effect group term estimate std.error statistic
fixed NA (Intercept) 1771.9066 41.84011 42.34947
ran_pars year sd__(Intercept) 118.0731 NA NA
ran_pars Residual sd__Observation 817.8890 NA NA 

# les variations des APC par année
tidy(model1_fit, effects = "ran_vals") %>%  # on calcule les valeurs des effets aléatoires
    mutate(estimate = 1772 + estimate,  # pour rajouter l'ordonnée à l'origine pour l'ensemble de la population
           ymin = estimate - 2 * `std.error`, # pour calculer des intervalles à 95%
           ymax = estimate + 2 * `std.error`) %>% 
    arrange(estimate) %>% 
    mutate(level = as_factor(level)) %>% 
    ggplot(aes(x = level, y = estimate)) +
    geom_pointrange(aes(ymin = ymin, ymax = ymax)) +
    geom_hline(yintercept = 1772, linetype = 2) +
    coord_flip()


  • Deuxième modèle : random intercepts par année et trend dépend de la discipline

Par rapport au premier modèle, on ajoute la discipline comme variable explicative du modèle - qui impacte donc l’angle de la pente de régression (mais celle-ci est la même par cluster, c’est-à-dire par année).

## Tentons un 2e modèle, dans lequel on suppose que le montant des APC dépend de la discipline
model2_fit <- model1_spec %>% 
    fit(amount_apc_EUR_trusted ~ bso_classification + (1 | year), data = train) # l'effet de la discipline sur l'APC est supposé le même quelle que soit l'année

# vue générale du modèle
summary <- tidy(model2_fit)
kable(summary, format = "html") %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), 
                full_width = F, 
                html_font = "sans-serif")
effect group term estimate std.error statistic
fixed NA (Intercept) 1849.35886 43.361568 42.649723
fixed NA bso_classificationChemistry -197.16272 15.207234 -12.965061
fixed NA bso_classificationComputer and information sciences -553.73092 19.908469 -27.813838
fixed NA bso_classificationEarth, Ecology, Energy and applied biology -232.68578 9.258577 -25.131917
fixed NA bso_classificationEngineering -486.09551 27.113104 -17.928434
fixed NA bso_classificationHumanities -882.32097 34.474890 -25.593148
fixed NA bso_classificationMathematics -925.75223 34.972179 -26.471106
fixed NA bso_classificationMedical research -13.18786 6.038334 -2.184023
fixed NA bso_classificationPhysical sciences, Astronomy -227.31944 10.420008 -21.815669
fixed NA bso_classificationSocial sciences -417.68790 34.692059 -12.039871
ran_pars year sd__(Intercept) 122.01518 NA NA
ran_pars Residual sd__Observation 804.11607 NA NA 


# les variations des APC par année - ici les valeurs sont calculées pour la discipline "Social"
tidy(model2_fit, effects = "ran_vals") %>%  # on calcule les valeurs des effets aléatoires
    mutate(estimate = 1852 + estimate,  # pour rajouter l'ordonnée à l'origine pour l'ensemble de la population
           ymin = estimate - 2 * `std.error`, # pour calculer des intervalles à 95%
           ymax = estimate + 2 * `std.error`) %>% 
    arrange(estimate) %>% 
    mutate(level = as_factor(level)) %>% 
    ggplot(aes(x = level, y = estimate)) +
    geom_pointrange(aes(ymin = ymin, ymax = ymax)) +
    geom_hline(yintercept = 1852, linetype = 2) +
    coord_flip()


  • Troisième modèle : random intercepts par année et trend dépend de la discipline par année

Par rapport au deuxième modèle, on autorise la pente liée à la discipline, à varier selon l’année (more random effects).

### Hypothèses

# Test d'homogénéité de la variance (H0)
bartlett.test(amount_apc_EUR_trusted ~ interaction(year, is_covered_by_couperin), data = bso_pop)
Bartlett test of homogeneity of variances

data: amount_apc_EUR_trusted by interaction(year, is_covered_by_couperin) Bartlett’s K-squared = 1959.9, df = 15, p-value < 2.2e-16

    #-->use weights arguments from lme()

# Test d'autocorrélation des résidus

# Taille des sous-populations
table(bso_pop$is_covered_by_couperin) # échantillons déséquilibrés

FALSE TRUE 102667 16947

table(bso_pop$year) # de 8184 à 26054 (peu d'obs pour 2013-2014 car suppression des NA)

2013 2014 2015 2016 2017 2018 2019 2020 8184 9537 11251 12978 14839 16943 19828 26054

    # 16 clusters
    #--> passer en bayésien


III - Imputation des valeurs manquantes de Y : 2013-2020

Homogénéité des articles avec Y inconnus

# Filtre sur les Y is_french_CA manquants
sample <- data %>% filter(is.na(is_french_CA_bso_wos_openalex_single_lang))


IV - Prévisions Y : 2021-2030

 

Document sous licence ouverte